home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / lin_alg.lha / lin_alg / matrix2.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  4KB  |  118 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Linear Algebra Package
  6.  *
  7.  *        Basic linear algebra operations, level 2
  8.  *           Matrix transformations and multiplications
  9.  *               of various types
  10.  *
  11.  ************************************************************************
  12.  */
  13.  
  14. #include "LinAlg.h"
  15. #include <math.h>
  16.  
  17. #pragma implementation
  18.  
  19. /*
  20.  *------------------------------------------------------------------------
  21.  *                Transpositions
  22.  */
  23.                 // Transpose a matrix
  24. Matrix Matrix::transpose(void) const
  25.      return tm(ncols+col_lwb-1,nrows+row_lwb-1,col_lwb,row_lwb)
  26. {
  27.   is_valid();
  28.   register REAL * rsp = elements;    // Row source pointer
  29.   register REAL * tp = tm.elements;
  30.  
  31.                 // Matrix tm is traversed in the natural
  32.                   // column-wise way, whilst the source matrix
  33.                 // is scanned row-by-row
  34.   while( tp < tm.elements + tm.nelems )
  35.   {
  36.     register REAL * sp = rsp++;        // sp = @ms[j,i] for i=0
  37.     while( sp < elements + nelems )    // Move tp to the next elem in the col,
  38.        *tp++ = *sp, sp += nrows;    // sp to the next elem in the curr row
  39.   }
  40.   assert( tp == tm.elements + tm.nelems && rsp == elements + nrows );
  41.  
  42. }
  43.  
  44. /*
  45.  *------------------------------------------------------------------------
  46.  *            General matrix multiplications
  47.  */
  48.  
  49.                 // Compute target = target * source
  50.                 // inplace,
  51.                 // target being the matrix the operation is
  52.                 // applied to.
  53. Matrix& Matrix::operator *= (const Matrix& source)
  54. {
  55.   is_valid();
  56.   source.is_valid();
  57.  
  58.   if( row_lwb != source.col_lwb || ncols != source.nrows ||
  59.       col_lwb != source.col_lwb || ncols != source.ncols )
  60.     info(), source.info(),
  61.     _error("matrices above are unsuitable for the inplace multiplication");
  62.  
  63.   REAL one_row[ncols];            // One row of the old_target matrix
  64.  
  65.   register REAL * trp = elements;    // Pointer to the i-th row
  66.   for(; trp < &elements[nrows]; trp++)    // Go row-by-row in the target
  67.   {
  68.     register REAL *wrp, *orp;               // work row pointers
  69.     for(wrp=trp,orp=&one_row[0]; orp < &one_row[ncols];)
  70.       *orp++ += *wrp, wrp += nrows;        // Copy a row of old_target
  71.  
  72.     register REAL *scp=source.elements;        // source column pointer
  73.     for(wrp=trp; wrp < elements+nelems; wrp += nrows)
  74.     {
  75.       register double sum = 0;            // Multiply a row of old_target
  76.       for(orp=one_row; orp < &one_row[ncols];)    // by each col of source
  77.     sum += *orp++ * *scp++;            // to get a row of new_target
  78.       *wrp = sum;
  79.     }
  80.   }
  81.  
  82.   return *this;
  83. }
  84.  
  85.  
  86.                 // Create matrix C such that 
  87.                 // C = A * B
  88. Matrix operator * (const Matrix& A, const Matrix& B)
  89.      return C(A.nrows,B.ncols,A.row_lwb,B.col_lwb)
  90. {
  91.   A.is_valid();
  92.   B.is_valid();
  93.  
  94.   if( A.ncols != B.nrows || A.row_lwb != B.col_lwb )
  95.     A.info(), B.info(),
  96.     _error("matrices above cannot be multiplied");
  97.  
  98.   register REAL * arp;            // Pointer to the i-th row of A
  99.   register REAL * bcp = B.elements;    // Pointer to the j-th col of B
  100.   register REAL * cp = C.elements;    // C is to be traversed in the natural
  101.   while( cp < C.elements + C.nelems )    // order, col-after-col
  102.   {
  103.     for(arp = A.elements; arp < A.elements + A.nrows; )
  104.     {
  105.       register double cij = 0;            
  106.       while( arp < A.elements + A.nelems )    // Scan the i-th row of A and
  107.     cij += *bcp++ * *arp, arp += A.nrows;    // the j-th col of B
  108.       *cp++ = cij;
  109.       bcp -= B.nrows;                // Return bcp to the j-th col
  110.       arp -= A.nelems - 1;            // arp points to (i+1)-th row
  111.     }
  112.     bcp += B.nrows;            // We're done with j-th col of both
  113.   }                    // B and C. Set bcp to the (j+1)-th col
  114.  
  115.   assert( cp == C.elements + C.nelems && bcp == B.elements + B.nelems );
  116. }
  117.  
  118.